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Abstract 

Comparison between three different numerical techniques for solving a cou- 
pled channel Schrodinger equation is presented. The benchmark equation, which 
describes the collision between two ultracold atoms, consists of two channels, 
each containing the same diagonal Lennard- Jones potential, one of positive and 
the other of negative energy. The coupling potential is of an exponential form. 
The methods are i) a recently developed spectral type integral equation method 
based on Chebyshev expansions, ii) a finite clement expansion, and iii) a combi- 
nation of an improved Numerov finite difference method and a Gordon method. 
The computing time and the accuracy of the resulting phase shift is found to be 
comparable for methods i) and ii) , achieving an accuracy of ten significant fig- 
ures with a double precision calculation. Method iii) achieves seven significant 
figures. The scattering length and effective range are also obtained. 
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I. Introduction 



The collision between two ground state atoms at low (/zK) temperature poses 
challenging computational problems which can be summarized by the words 
" long range" and " coupled channels" . The former problem arises from the fact 
that in the collision process the atomic clouds "polarize" each other, leading 
to long ranged dispersion potentials.il! The latter problem arises from the fact 
that the internal hyperfine structure of the atoms leads to a set of coupled 
Schrodinger equations that describe the transitions which an atom in the inci- 
dent channel can make to many of these hyperfine states. Accurate calculations 
of the scattering properties and wave functions of the atom-atom collision are 
crucial, since the macroscopic shape of the Bosej-Einstein condensate as well 
as the lineshapes of the photoassociation spectraH depend sensitively on these 
quantities . 

The lower the incident energy of the atoms, the larger are the distances for 
which the potentials affect the phase shifts. This can be seen from the WKB 
approximation to the phase shift which contains integrals over the local wave 
length k(x) = [2fj,(E - V(x))} 1/2 of the form 



f k{x)dx ~ f kodx — f V(x)dx, 

J J 2feo J 



since the first significant term in the expansion above contains the ratio of 
the potential to the asymptotic wave number. Further, when negative energy 
channels are coupled to the positive energy incident channel, it can become 
difficult to enforce the appropriate decaying wave function boundary condition, 
with resulting loss of stability, depending on the algorithm used. 

There are various calculational methods available for dealing with this scat- n 
tering problem: modified Numerov, Gordon'scl,cl a finite element method (FEM),Ea 
and a recently developed method that consist of replacing, the coupled differ- 
ential equations by equivalent integral equations_(IEM)j3.0 In addition, there 
are more sophisticated finite difference methodsa; we, however, will not con- 
sider these methods here since they are not as widely in use. Methods in- 
volving the representation pJ a continuous function by a finite set of sampling 
points have been discussedJa One such method led to the mapped Fourier grid 
methodJi3 which has been employed for the calculation of the collision between 
cold atoms. El The interaction of cold atoms with surfaces has also been dis- 
cussed J13 including how to implement boundary conditions.. 

Depending on the degree of accuracy required and the ease of performing the 
calculation, any one of these methods may be the most suitable for a particular 
situation. It is nevertheless of interest to compare these methods with each other 
as far as accuracy, stability and numerical complexity are concerned. It is the 
purpose of this paper to compare the three methods, by numerically evaluating 
a benchmark test case described below. 
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"Numerical computational stability" is a many- faceted concept. It mani- 
fests itself through the degree of accuracy obtained. There are at least three 
ingredients: 

• i) How the numerical truncation error of the algorithm is offset by the 
presence of the accumulation of round-off errors. The larger the number 
of mesh points the smaller the truncation error, but the larger is the cor- 
responding overall round-off error. Different algorithms strike a different 
balance between these two errors. 

• ii) The sensitivity of the final result to the errors in the input data, such as 
the potentials, masses, etc. This sensitivity is expresses by the "condition" 
of the formulation model. 

• hi) How the asymptotic boundary conditions are achieved, both in the 
open and the closed channels. For the latter, the growing solutions that 
contaminate the decaying solutions have to be eliminated. Each algorithm 
proceeds by a different method. For example, in the IEM we introduce 
scaling factors in each partition that prevent the unwanted solution. In 
the finite difference methods, one has to integrate inward and outward and 
then match at some intermediary distance. Here the matching matrix can 
introduce errors. With the finite elements, the closed channel solutions are 
forced to be zero at the final matching point, automatically eliminating 
the unwanted growing solutions. 

It is not the purpose of this paper to investigate in detail the stability prop- 
erties of the three algorithms described in the present study. For that, a detailed 
comparison of the numerical solution with an exact solution for an artificial test 
case would have to be performed. Rather, we here attempt to obtain some nu- 
merical evidence for the degree of stability of the three methods for a realistic 
example. 



II. The Test Case 

The test case we have chosen consists of a model calculation that captures 
the essence of the collision of two ultra-cold 2 S alkali atoms. The character- 
istic feature of such collisions is that in going from large to small internuclear 
separations a change of coupling schemes occurs. Asymptotically, the hypcrfinc 
structure of the individual alkali atoms dominates, while at short internuclear 
separations the molecular X 1 T lg and o 3 E„ potentials dominate. A basic under- 
standing of the properties of these two Born-Oppcnhcimer potential curves can 
be found in any text book on quantum mechanics that discusses the electronic 
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structure of a H2 molecule. Here it is sufficient to realize that they have an iden- 
tical long-range attractive van der Waals behavior and are split exponentially 
via an exchange mechanism at shorter internuclear separations. 

The simplest multi-channel potential that captures this physics is thus a two 
channel model. Taking into account that for small collision energies the nuclear 
rotation can be safely ignored, the Hamiltonian is conveniently parametrized as 



where the reduced mass /i = M/2, E^f is the asymptotic splitting between the 
two channels, and £ is the total energy in the system. The r-dependent poten- 
tials of our test problem are the Lennard- Jones potential Vuir) — Cu/r 12 ~ 
Ce/r 6 and an off-diagonal exchange coupling given by Ae~ br . The functions ipp 
and ipN describe the wavefunction for the open and closed channel, respectively. 
Notice that the zero of energy is located at the lowest of the two asymptotes. 

Obviously this Hamiltonian is set up in terms of the atomic basis. At large in- 
ternuclear separation the Hamiltonian reduces to a diagonal matrix. In fact, the 
Hamiltonian would be diagonal for all internuclear separation if the exponential 
off-diagonal potential were absent. It turns out that for internuclear separations 
where this term is large compared to £hf it is informative to calculate the adia- 
batic potentials by diagonalizing the potential term of the Hamiltonian at each 
internuclear separation. At shorter distances the resulting potentials correspond 
to a very deep X 1 E S and a shallow a£ u potential, to a good approximation. 

For two colliding ultra-cold 2 S Na atoms—realistic values of the constants are 
M = 22.9897680 amu, C 6 =1472 a.u.(a ) 6 H C*i 2 =38 x 10 6 a.u.(a ) 12 , A=2.9 
a.u., b=0. 81173 ao, and Em— 0.2693-10 -6 a.u. This choice of Em is approx- 
imately equal to the atomic hyperfine splitting of the 2 S Na atom. The total 
energy E = 3.1668293 x 10~ 12 a.u. corresponds to a temperature of 1 /zK. Since 
E <C Em j the energy in the second channel is negative, i.e., only one of the two 
channels is asymptotically accessible. In the above a.u. stands for atomic units, 
and ao is the Bohr radius. 

The conversion into entirely ao units is achieved by dividing the above equa- 
tion by 2/i/h 2 . One obtains 



where r is in units of ao and the potential and energy matrices, V and £ re- 
spectively, are in units of (ao)~ 2 . The conversion of a quantity in a.u. units 
to (ao) -2 units is achieved by multiplying the former by /i = 22.989768 x 
1822.888506(a.u.) -1 (a ) -2 . The potential matrix is 




(1) 





V = 



V u 
U V 



(3) 
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where V = V X fj,, U = U X (a, and the energy matrix is 



Here k and k, the wave numbers in each channel, are given by k = \JE x fi and 
k = y/Ehf x fi — k 2 . In our numerical example, the corresponding values are 
k = 3.643004224146145 x 10- 4 (a )" 1 and k = 0.1062338621818394(a )- 1 . The 
wave function is normalized so that asymptotically it becomes 

( \ ~ ( Sm (^ r ) + Ki cos(fcr) 
\ ipN J ~ V K 2 exp(-nr) 

where K\ and Ki are two elements related to the real scattering R matrix, in 
terms of which the phase shifts can be obtained. 

In this model, the diagonal potential extends to considerably larger distances 
than the coupling potential. Further, between 6 and 10 ao the diagonal potential 
is very deep leading to many oscillations in the wave functions. For example, 
near 5.5ao the local wave length A in both channels is ~ 0.25ao, near 8.5ao 
A ~ 1.2ao, with smaller ripples superimposed, and near 20ao the local wave 
length has increased to ~ 4ao At distances less than 4ao the repulsive portion of 
the potential becomes very large making the wave function very small. In order 
to allow for the singularity of the diagonal potentials near the origin, a parameter 
R cu t is defined, and the wave functions are set to zero in the interval [0, R cu t\- A 
value of R cu t = 4.0ao is found to be satisfactory. In addition, the calculation is 
carried out to a maximum radius, i? maX j beyond which all potentials are set equal 
to zero. In our calculations, R max is set equal to 500a o . When R max is increased 
further, the values of K\ and K2 still change beyond the 6 th. significant figure, 
as can be seen from the Table in Appendix, even though the Lennard-Jones 
potential is less than —3.95 x 10 -9 a^ 2 . The large effect on the phase-shift 
produced by such a small potential is due to the occurrence of the factor 1/fc ~ 
2.75 x 10 3 in the integrals involving the potential tail, as was already pointed 
out in the introduction. Rather than numerically calculating such changes, it is 
preferable to employ perturbation methods, which are described in Appendix 2. 
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A. The Integral Equation Method. 

In this method the differential equation + ip — Vip is transformed 
into the Lippmann-Schwinger integral equation 

V>(r) =F(r)+ / g {r,r')V{r')dr', 



where F(r) is a undisorted wave function, like (sin(fcr),0) , and Go(r, r') is the 
undistorted Green's function matrix.Q A motivation for such an approach is 
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that the solutions of integral equations have better numerical stability than the 
solutions of differential equations. One common objection to the use of inte- 
gral equations has been that the solution leads to large matrices which are not 
sparse and hence require substantially larger amounts of computing time than 
the sparse matrices of differential equations. This objection was overcome in our 
integral equation method (IEM) by dividing the whole radial interval into par- 
titions. The integral equations in each partition lead to dense matrices of small 
dimension, but the matrix that combines the local solutions into the global one, 
albeit of large dimension, is sparse. It should be noted that this latter property 
is valid only in configuration space, because only in this space do the Green's 
functions have the required semi-separable nature. In the present version of the 
IEM method the (variable) size of each of the partitions is determined in terms 
of two parameters NL and e as follows. In each radial region a local wave length 

in channels 1 and 2 is obtained as 2tt/ \J\k 2 — V{r)\, and 27r/y| — k 2 — V(r)\. 
The smaller of the two local wavelengths is taken, and the size of the partition 
in that region is determined such that there are a given total number NL of 
Chebyshev points per local wave length. Allowing for the fact that in each 
partition there are 16 Chebyshev points, the average length of a partition for 
a given local wave length A is A x 16/NL. The length of each partition is sub- 
sequently readjusted using the tolerance parameter e as follows. According to 
the IEM method,Q in each partition two sets of "local" functions are calculated 
in terms of which the global function ?/> is obtained as a linear combination. 
The accuracy of each of the local functions can be determined by the size of the 
coefficients of the highest order Chebyshev polynomials. If the relative accuracy 
of the local functions in a given partition is larger than e. then that partition is 
divided in half, and the testing is continued. If the initially chosen value of NL 
is too small, then the initial partitions are too large, and many of the partitions 
are subsequently reduced by the e criterion. In this case the final number of 
partitions M becomes larger than their initial value. If the chosen value of NL 
is too large, then most of the partitions are unnecessarily small, and the value of 
M is too large, leading to a larger accumulation of roundoff errors for the final 
elements of the i-L -matrix. An exception is the interval [0, R cu t]- This interval 
is considered as one partition, containing a total of 16 Chebyshev points. This 
is sufficient since the wave function is very small in this region (less than the 
desired accuracy for the values of K ), and the values of K were found to be 
stable to 11 significant figures as R cu t was varied below 4.0ao . 

In summary, for a given value of e, the value of NL was varied until the 
smallest number of partitions M was obtained. An example is given in the 
table below. 



Table 1. Values of K\ and number of partitions M as a 
function of NL for the tolerance e. = 1CP 9 
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NL 


M 


Ki 


K 2 


10 


150 


-0.31233398338809 


6.5761303971514 


20 


153 


-0.31233398339572 


6.5761303973071 


30 


144 


-0.31233398339229 


6.5761303972290 


40 


154 


-0.31233398339070 


6.5761303972039 


50 


177 


-0.31233398338870 


6.5761303971639 



From this table one can find a value of K\ — —0.3123339834 and K 2 = 
6.576130397 which are stable to ten significant figures. For values of the tol- 
erance e between 10~ 13 and 10~ 3 a good compromise value of 10 for NL was 
found. The corresponding values of M and the corresponding accuracy of K\ 
are listed in Table 2 for several values of e. 



Table 2. Accuracy for K\ and number of partitions M for a given 
value of the Tolerance e, with NL = 10 



e 


M 


K x 


# of Sign. Figs. 


IQ-' A 


28 


-O.31243337402099 


3 


10~ 3 


30 


-0.31233467247746 


5 


10~ 5 


67 


-0.31233398457315 


8 


io- 7 


106 


-O.31233398370637 


8 


io- y 


158 


-0.31233398338809 


10 


io- n 


214 


-0.31233398338869 


10 


io- 13 


574 


-0.31233398338534 


10 



From this table it appears that beyond e = 10~ 9 the accumulation of roundoff 
errors begins to dominate, and -0.31233398339 is the best value of K\. The 
distribution of partitions for three tolerance parameters is shown in Fig. l.The 
increasingly large spacing of the partitions at the large distances is clear from 
the figure. In the vicinity of R ~ 50 the density of partitions is high because 
the negative energy channel has a turning point there. This shows that the 
local wave-length criterion alone would have been insufficient to determine the 
partition size. 

1. Scattering Length and Effective Range. 

As a further test of the stability of the IEM method, the Scattering Length 
a and Effective Range r e are investigated. They are obtained in the limit of 
small wave number k from the expression 

k/K x =-i+r e fc 2 + (3(fc 3 ). (6) 
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Figure 1: Various IEM partition distributions, described in Table 2. The y-axis 
represents the partition number i, and the x-axis shows the lower boundary 
of partition i. The more points in a particular radial interval, the smaller are 
the lenghts of the partitions in that interval. All three partition distributions 
started with the same number of mesh points per local wavelength (NL = 20), 
but were subsequently modified by the e-accuracy criterion, with e = 10~ 9 , 10~ 6 , 
and 10 -3 , respectively. The numbers above each curve represent the number 
of accurate significant figures achieved for the asymptotic constant Ki for each 
value of e. The large concentration of partitions in the vicinity of 50ao reflects the 
occurrence of a turning point in the negative energy channel near that distance, 
where accuracy would have been lost had the initial partition distribution been 
used. 
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The left hand side of the above equation is usually written as fccot^o, which is 
equal to k/K\. For a given choice of the truncation radius i? ma x the value of K\ 
is calculated for two different and small values of k and the values of a and r e 
are obtained from the two values of k/K\ in Eq. (||). 

However, as mentioned in the introduction, the values of the scattering ma- 
trix K depend on the choice of the truncation radius, the more so the smaller 
the value of k , because of the increasingly non-negligible contributions of the 
potential beyond i? m ax- Hence, in order to obtain a reliable value of the scatter- 
ing length, it is advantageous to first correct K\ for the contributions beyond 
i? max - This can be done quite simply by using first order perturbation theory, 
as is detailed in Appendix 1, and as will be demonstrated below. 

An example of the variation of the if's with i? max , obtained by the IEM, is 
given in the table below. All results in this section were obtained with NL = 
10, and e = 1CT 9 . . 



Table 3: Dependence of K\ and K 2 on R, 



^l(i?max) 


Xi(oo) 


K 2 (R max ) 




-.313705209 


-.312322902025 


6.62005410 


250 


-.312333983 


-.312323344009 


6.57613040 


500 


-.312324588 


-.312323343934 


6.57555968 


1000 


-.312324073 


-.312323343936 


6.57558157 


1500 


-.312323719 


-.312323343936 


6.57558741 


2000 



One can see from the table that the value of Ki(R max ) becomes mono- 
tonically less negative as i? ma x increases, while K 2 first decreases, and then 
increases for i? m ax > 1000. This behavior can be reproduced numerically by 
means of perturbation theory, described in the appendix. The usefulness of the 
perturbative correction is also demonstrated by the stability of the column de- 
noted as -fTi(oo), which contains the corrected values of K 1 (R ma ^,). The table 
shows that the perturbative correction increases the stability of K\ from 5 or 6 
significant figures to 11, yielding ^(00) = —.31232334394 

The stability of the values of a and r e will be described next. The two 
smallest values of k to be used in Eq. (||) were approximately 0.1152 x 10~ 4 and 
0.3643 x 10~ 6 , for which the third order term is smaller than the accuracy of 
the present method. Their values with and without the perturbative correction 
of Ki are shown in the table below. 



Table 4: Dependence of a and r e on R ! 
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a(-Rmax) 


a(oo) 


Ve (^max) 


r e (oo) 


^max 


852.0123407 


851.9817159134 


55.08319944 


55.1051720100 


500 


851.9849554 


851.9817157297 


55.10378197 


55.1051694472 


1000 


851.9837969 


851.9817157354 


55.10138602 


55.1051696827 


1500 


851.9829235 


851.9817157362 


55.10046968 


55.1051693964 


2000 


851.9824507 


851.9817157362 


55.10031183 


55.1051691577 


2500 



The table shows that the stability of the scattering length is increased from 
5 to 11 significant figures by the perturbative correction, yielding a(oo) = 
851.98171574. Likewise, the stability of the effective range is increased from 
5 significant figures to 8, yielding r e (oo) = 55.105169. 

B. The Finite Element method 

The second method employs the non-iterative eigenchapael0 variant of the 
R-matrix method first introduced by Wigner and Eisenbud.ll3 The eigenchannel 
R-matrix method solves the Schrodinger equation within a finite reaction vol- 
ume n of configuration space, subject to constant normal logarithmic derivative 
boundary conditions on the surface E of f2. The collisional properties of the 
system, typically represented in terms of an S'-matrix, are easily obtained once 
the normal logarithmic derivative b — [d^ /vitialn)^^ 1 is calculated. 

One can obtain the following eigenvalue equation for the normal logarithmic 
derivative b on the surface 5x2 

Tc = bAc (7) 

We solve this equation using the finite-element method (FEM).I The FEM di- 
vides the radial domain into N sectors (or partitions) and within each sector 
defines a local basis in much the same spirit as the IEM above. The local ba- 
sis functions, however, are fifth-order Hermite interpolating polynomials rather 
than Chebyshev polynomials. The Hermite interpolating polynomials Uk(x n ), 
fc=l-6, are non-zero only in sector n. Here, rescalcd variable defined on 

the interval [—1,1] that is related to the physical internuclear separation R in 
sector n through an appropriate linear transformation. The six basis functions 
are defined through the following boundary conditions: 

Mfc(-l) = $ik «fc(0) = S 3k u fc (l) = S 5k , . 

u' k {-l) = 5 2k u' k (0)=6 4k u' k (l) = 6 6k ■ W 

A multi-component radial wave function can then be represented by the follow- 
ing expansion on the piecewise polynomials Ui(x): 

(9) 

i = { m , /c , n } 
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The set i contains the basis function index k, the channel index m, and the 
sector index n. The coefficients Ci are to be determined subject to a continu- 
ity constraint on each channel component and its first derivative across sector 
boundaries. The assembly stage of the calculation is thus performed prior to 
the main calculation, in contrast to the IEM where it is performed after. In ad- 
dition, channel boundary conditions can be imposed quite simply by setting the 
value of the appropriate coefficient. For instance, the closed channel function is 
forced to be zero on the surface by setting the coefficient of % in the last sector 
to zero. 

In the finite element representation, the matrix elements of T and A in Eq.(Q) 
are given by 

Tij = 2n J Ui(x n )(E- H)Uj(x n )a n dx n -S m , m >S niN 5 kt5 8 k ^ 6 /a n (10) 

and 

Ajj = <5m,m"5n,iV^fc,5<5/c',5 , (H) 

respectively, where a n = (i? n +i — -Rn)/2. Here i = {m, k, n} and j = {m', fc', n}, 
and H is the Hamiltonian of the system. Because the basis functions are non- 
zero only within a given sector, the corresponding matrices have roughly a 
block diagonal structure. Each sector has an associated block which is coupled 
only to its nearest neighbors through the continuity conditions. Note that the 
overlap matrix A, whose elements are given by surface integrals over the basis 
functionsJi-3 is particularly simple in the FEM representation. It contains only 
m a non-zero elements (all equal to one), where m is simply the number of 
channels open (i.e., E > V m ) on the surface S. The T matrix is symmetric and, 
as mentioned, can be constructed in a banded format. The integrals representing 
the matrix elements of T_ are also particularly simple in the FEM representation. 
In fact, except for the integral over the interaction potential, all integrals can 
be evaluated analytically once and for all before hand, significantly decreasing 
the CPU time necessary to construct the matrix. 

At this stage, we are left with a banded generalized eigensystem to solve that 
typically has large dimensions but has at most only m non-zero eigenvalues. 
Although there are standard linear algebra packages which could solve these 
equations directly, implementing an efficient, general method would be difficult 
since the non-zero eigenvalues can range— between — oo and +00 and all mo 
of them are needed. It has been shown,li-3 however, that by partitioning the 
matrices according to whether the basis functions are non-zero (open = o) or 
zero (closed = c) on E, i.e. by writing 



r cc r co \ f c c \ \ / c c 

r oc r°° H c° j I A°° ) I c° 



(12) 



Eq. (0) can be reduced to a small (m c x m ) eigensystem 

£l 00 c° = bA 00 c° (13) 
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where Q 00 = T°° - r oc (r cc )- 1 r c0 . The main computational burden is thus 
shifted to constructing £2°°. Since T_ cc is a large banded matrix in the FEM 
representation, this is most efficiently accomplished by solving the set of linear 
equations 

T CC X = T co . (14) 

The solution X = (r cc )-ir co thus provides the needed matrix inverse. T cc 
has dimensions AMNxAMN, but with a half-bandwidth (number of non-zero 
diagonals above the main diagonal) of only 6M — 1 resulting from the FEM 
representation. r co has dimensions 4MNxm . We use standard LAPACKO 
routines to solve equations [l^ and [l4|. The eigenvalues b and eigenvectors c° 
completely specify the linearly independent solutions of the Schrodinger equa- 
tion on the surface S. This information is generally packaged in terms of a 
i?-matrix 

Rmm' = ^ Z mpbp 1 Z f3 ^ l , (15) 

where the columns of Z are given by the eigenvectors c° . The ^-matrix is then 
obtained through simple matrix manipulations involving only the i?-matrix and 
the two linearly independent solutions of the asymptotic form of the Schrodinger 
equation.li.3 

For the two-channel test problem described above, Eq. ( [l4| ) is a 8Nx8N 
matrix equation with a half-bandwidth of 11, and only one solution is required 
since there is only a single open channel, m —l. The present formulation of the 
R-matrix method focuses the computational effort on finding only the relevant 
scattering information. Thus, only the coefficient in the open channel K\ is 
obtained. The results are shown in the table below as a function of the number 
of sectors N used. 



Table 5. Results for the Finite Element Method 



Ki 


N 


-0.312345739663914 


400 


-0.312334009008278 


800 


-0.312334008856244 


1600 


-0.312334008759115 


3200 


-0.312334006921697 


6400 



C. The Gordon Algorithm. 

The Gordon algorithm^ is a well-established numerical method to solve 
for the scattering solutions of a set of N c coupled radial Schrodinger equations. 
Similar to the FEM discussed in the previous section, this algorithm is local in 
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the sense that the wavefunction is propagated from R to R + dR using only the 
wavefunction at R. The version of the Gordon method used here is described 
in Ref..O The main difference from the original Gordon method is in replacing 
the Airy functions, which correspond to the reference solutions for linearized 
potentials, by singr or cos qr in each interval dR. Here the quantities q are the 
local wavenumbers of each channel in the interval dR . They are obtained by 
replacing the potential matrix in interval dR by an averaged constant potential, 
and then diagonalizing this potential. If the local energy is negative then exp(qr) 
or exp(— qr) are used instead. The wave function in the diagonalized interval 
dR is represented in each channel n by \t n = A n (r) sing„r + B n {r) cosg n r,and 
the vector of the coefficients A and B are calculated by solving a system of first 
order coupled differential equations involving the difference between the true 
and the averaged potential matrix in interval dR. In the intervals where some 
of the local energies are negative (i.e., some of the q's are imaginary), the unde- 
sired exponentially increasing functions are minimized by a " triangularization" 
method developed by Gordon. Analytic connection formulas between R and 
R + dR determine A and B, and thus the solution is propagated towards the 
final matching point. Consequently no large linear system needs to be evaluated 
in the Gordon algorithm and hence the method is not memory limited. 

The physical boundary conditions for the solution are obtained by first calcu- 
lating N c linear independent solutions from different initial conditions near the 
origin, and then constructing appropriate linear combination of these "math- 
ematical" solutions. The stepsize dR is calculated on the fly. After a fixed 
number of steps (say ten) the percentage change of the A and B coefficients 
between intervals n and n + i, where i ~ 10, is compared to an accuracy cri- 
terium. If the variation is too large, the step size is repeatedly halved until the 
accuracy criterion is met. On the other hand if the coefficients hardly changed 
the step size is doubled. Consequently, if the potential is well approximated by a 
constant, large steps are taken. Closed channels are removed from the propaga- 
tion when the corresponding amplitude of the channel wavefunction has become 
smaller than a threshold parameter. The influence of these components on the 
scattering properties of the asymptotically accessible channel is then negligible. 

Our results with this algorithm for R cut =A ao and R max =5Q0 a are sum- 
marized in Table 6. The columns describe the fractional change of the A and B 
coefficients, the number of steps, and K\, respectively. It is immediately clear 
that the Gordon method uses many more steps than the other two numerical 
methods. As discussed before this is not crucial as it is not necessary to store 
the wavefunction at every step in order to propagate the wavefunction. 

Table 6. Results for the Gordon Method 
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Fract'l 


No of 




Change 


steps 


K x 


0.1000 


7025 


-0.312348088 


0.0500 


7632 


-0.312340564 


0.0250 


8353 


-0.312337479 


0.0100 


9796 


-0.312335340 


0.0050 


11402 


-0.312334854 


0.0020 


14454 


-0.312334070 


0.0010 


18179 


-0.312334067 


0.0005 


25502 


-0.312334020 



III. Discussion 

A comparison between the three methods is displayed in Table 7, which lists 
the number of mesh points required to achieve a certain accuracy for K\ . 



Table 7. Comparison between three methods. 







No of 


No of 


Total No 




K x 


partit's 


Pts/part. 


of Points.. 


IEM 


-0.3123339834 


153 


16 


2,448 


FEM 


-0.312334009 


800 


4 


3,200 


Gordon 


-0.3123340 


25,502 


1 


25,502 



The finite element (FEM) and the integral equation (IEM) methods are 
nearly identical in performance. Both can easily adjust the size of the partitions 
to the local conditions of the potentials; both give stability of at least nine 
significant figures, and both use approximately the same number of mesh points. 
Their numerical complexity is also comparable since a CPU-time test shows that 
both use approximately the same computer time for the case tested here. The 
IEM and the FEM differ in the 8th significant figure (by 2.6 x 10~ 8 ) for K x . 
The reason for this difference is not known, but could be related to the fact 
that both the FEM and the Gordon methods have not yet fully converged as 
the number of partitions is increased, as can be seen from Tables 5 and 6, and 
in Figs. 2 and 3. These figures show the value of K\, from the fifth significant 
figure onward, as a function of the total number of mesh points for each of the 
three methods in progressively larger detail. The somewhat slower convergence 
for the FEM may be due to the fact that the algorithm for determining the 
partition size, especially in the vicinity of a turning point, is not as refined for 
the FEM as for the IEM. 
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Figure 2: Comparison of the rate of convergence of the numerical value of K\ as 
a function of the total number of mesh-points N. On the vertical axis is plotted 
the value of K\ from the fifth significant figure onwards. This is accomplished 
by subtracting -0.312300000 from each value of K\ and multiplying the result 
by 10 4 . The values of K x are taken from Tables 2, 5, and 6 for the IEM, FEM 
and Gordon methods. They are represented by solid circles, open circles and 
squares, respectively 
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Figure 3: Same as Fig. 2 with a larger magnification. 
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The Gordon method is simpler to implement than the two other ones and 
gives a respectable accuracy of seven significant figures, albeit at the expense 
of a much larger number of mesh points, which in turn leads to a larger ac- 
cumulation of roundoff errors. An improved Numerov method gave only four 
significant figures {K\= -0.31233(2)) and is not mentioned further. In the IEM, 
the boundary conditions are built in automatically via the Green's functions, 
while in the present form of the FEM the solutions with exponential growth are 
eliminated by forcing the closed channel component to be zero on the surface 
S. In a separate study of two coupled equations]!! it is shown that there are 
situations in which the conventional Numerov method has severe difficulty in 
obtaining the correct asymptotic boundary condition, while both the IEM and 
the FEM do not. 

It should be clear that the accuracy achieved in this benchmark calculation 
is not directly relevant for physical applications, since the potentials are usually 
known only to low accuracy. Instead, the accuracy achieved is to be construed 
as a measure of the stability of the method, which gives an indication of how 
well the method is expected to perform under more complex situations, such 
as when many channels are involved, or when the range of the interaction is 
excessively large, or at very high energies where many oscillations in the wave 
functions are present. 

In conclusion, all three methods performed well in providing a numerical 
solution to the coupled channel test case examined here. This case could serve 
as a benchmark calculation for testing additional methods, since the effective 
range and the scattering length are also calculated. Further comparison of 
various methods under more complex conditions would be desirable in order to 
determine the conditions under which a particular one of the methods would be 
preferable. 



Acknowledgments: E.T. would like to acknowledge useful discussions with 
Dr. Fred Mies, who has implemented the original version of the Gordon method 
used by the authors. 

IV. Appendix 1. Perturbative treatment of the long-range 
corrections. 

As the value of i? max (this is the truncation radius beyond which all poten- 
tials are set to zero) is increased, the values of K\ and K2, defined in Eq. (||) 
change slightly because of the non-zero value of the potential at large distances. 
In order to extrapolate the K values to i? max = 00, it is preferable to include the 
long-range tail of the potential perturbatively rather than numerically. Showing 
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how this can be done is the purpose of this appendix. It is of course also-possi- 
ble to include the long range tails numerically, as has been done in Ref.ll3 The 
variation of the if's with large distances, as given in Table 3, can be reproduced 
numerically by means of perturbation theory, as will now be shown. 

For this section we assume that the radial distance is sufficiently large so that 
the coupling potential between the two channels is negligible. Thus, only the 
effect of the diagonal potential in each channel i on Ki needs to be considered. 
We also continue to assume that the angular momentum number I is zero. Both 
assumptions can removed by a generalization of the present discussion. 

A. The positive energy channel. 

We will denote two consecutive values of i? m ax by T\ and T2., respectively, 
and assume that the numerical calculation extends out to T\. We denote by 
Vi(r) the diagonal potential for r < T% and set Vi(r) = for r > T\. The 
corresponding wave function is denoted by ipi which is equal to the wave function 
tpp known numerically. Because Vi(r) = for r > T\, 

tpl(r) = sin(fcr) + Ki{T x ) cos(fcr). (16) 

Similarly, V^(r) is equal to the diagonal potential for r < T2 and V^(r) = for 
r > T2. The corresponding wave function, to be calculated perturbatively, is 
denoted by ip2- It obeys the equation (d 2 /dr 2 — V2 + k 2 ) -02 = 0, which can be 
written 

(d 2 /dr 2 + k 2 ) V> 2 = V'4> 2 r > T u (17) 

where 

V'(r) = V 2 {r) forr>T t (18) 
V'(r) = jorr<T x . 

Thus, V is the perturbative potential which vanishes outside of the interval 
P~i, T2]. Using for the inverse of the operator (d 2 /dr 2 + k 2 ) the Green's function 
integral expression in terms of sin(fcr<) x cos(fcr>), one obtains the most general 
form for tp2 

^(r) = asin(fcr) + (3 cos(fcr) + sin(/cr)E c (r) + cos(kr)T, s (r) (19) 

where ^ 

E c (r)=-i f 2 cos{kr')V'(r')iP2(r')dr' (20) 
k J r 

and r 

E.(r) = -\ I sm(kr r )V'{r r )^2{r r )dr'. (21) 
k J Tl 
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The coefficients a and (3 are obtained by matching ip 2 to ipx at r = Ti, i.e., setting 
the two functions and their derivatives equal. Remembering that E s (Ti) = 0, 
that (1/fc) {d^2/dr) r=Ti = a Sl +a (-/? + X c (Ti)) , and that (1/fc) (#i/dr) r=Ti = 
ci — Ki(Ti)sx, where s\ = sin(fcTi) and c\ = cos(fcTi), one obtains 

a = l-E.(Ti), /3 = K 1 (T 1 ). (22) 

Inserting the above result into Eq. (|l9f), one obtains for r >T 2 the result 

V> 2 (r) = (l-S c (T 1 ))sin(fcr) + (K 1 (T 1 ) + S s (T 2 ))cos(fcr). (23) 

The first order perturbation result for the above expressions is obtained by 
replacing ip 2 by ipi m the integrals for S c and E s in Eqs. above. Denoting the 
corresponding integrals from T\ to T 2 by I c and I s , respectively 

I c = -j f 2 cos{kr')V'{r')ilJx{r')dr' (24) 

fc JTi 

/ s = -i [ 2 sm{kr')V'{r')i) X {r')dr', (25) 

fc JTi 

one finally obtains the first order result for iTi(T 2 ) 

i^ 1} (T 2 ) - A ^ ) +/s ~ tfi(T x ) + I s + K^Tx) x Z c . (26) 

t -/c 

This is our final expression. The numerical evaluation was carried out ini- 
tially by using MATHEMATICAtJto evaluate the integrals involving products 
of circular functions and 1/r 6 . For example, if T\ and T 2 are set equal to 500 
and 2000ao respectively, and using for Cq the value stated above, one finds 
I c = -8.443 x 10~ 5 and I s = -1.6105 x 10~ 5 . Using for K 2 (500) the value from 
Table 3, one obtains from Eq. ( p6| ) the result 

i^ 1} (2000) - K 2 {b00) = 1.0265 x 10 -5 , 
which compares very well with the numerical result from Table 3 
Jf 2 (2000) - K 2 (500) = 1.0264 x 10~ 5 . 

The results for the scattering length and the effective range, described in 
the text, were calculated by a FORTRAN code which expresses the integrals 
/ [sin(£) /t]dt and / [cos(t)/t]dt in terms of the functions Ci and Si. The latter 
were called from the IMSL scientific library, and generalization to larger powers 
of t in the denominator (like t e in the present case) were obtained recursively 
through integrations by part. 
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B. The negative energy channel. 

For the negative energy channel a similar perturbative procedure will now be 
described. In what follows ip\ and ij) 2 now refer to the negative energy channel 
but otherwise they have the same meaning. The function tpi coincides with the 
numerical solution (/V for r < Ti, while </'2 for Tj < r is given by 

Mr) = lG 2 (r) - -F 2 (r) [ ' G 2 V'^ 2 dr' - -G 2 {r) [ F 2 V'^ 2 dr'. (27) 

Here F 2 (r) and G 2 (r) are equal, respectively to exp(— nr) and sinh(Kr). By 
setting ip 2 equal to ip± at r = Ti one obtains 

7 = T^FT (Wi) - - T 2 G 2 (r')V'(/)V 2 (r')dr'F 2 (T 1 ) ] , 



G 2 (Ti) 



Inserting the above into Eq. (p7|) evaluated at r — T 2 , replacing i\i 2 in the inte- 
grals in the above expressions by 

Mr) K 2 (Tx)G 2 (r) 
and dividing the result by G 2 (T 2 ), one obtains a preliminary value for K 2 (T 2 ) 

TV (Til 

^ 2 (T 2 ) - Jf 2 (T x ) - ^rWe + he 



where 



h = r G 2 (r r )V'(r')G 2 (r')dr' 



H 



hs = Sll f 1 F 2 (r')V'(r')G 2 (r')dr' (28) 

The final value of K 2 (T 2 ) is obtained by dividing K 2 (T 2 ) by the same normaliz- 
ing factor (1 — I c ) which was required to normalize the wave function in channel 
1 so that the coefficient in front of sin(fcr) be equal to 1. The final expression is 

K 2 (T 2 ) = ^2jL ~ K2{Tl){l + Ic) _ g^/e + hs (29) 

Using for Ti and T 2 the values 500a and 2000ao respectively, numerical eval- 
uation of Eq. @ gives # 2 (T 2 ) - K 2 (Ti) ~ -5.445 x 1CT 4 while the numerical 
value obtained from Table 3 gives ~ —5.430 x 10~ 4 
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